#!/usr/bin/perl

#-------------------------------------------------------------------------------
# Script to get components from a spectClean.exe spectrum file
#-------------------------------------------------------------------------------


# Variables
#===============================================================================
$PI         = 3.141592654;
$spectrum   = "clean.dat";
$beamParamters = "beamParameters.dat";
$startFreq = 0.0;
$endFreq = 1E2;
#$endFreq = 1.0;
#===============================================================================



# Get and parse command line variables
$numArgs = $#ARGV + 1;
print STDERR "\n";
print STDERR "------------------------------------------------------\n";
print STDERR "ARG1 (DIR)    = $ARGV[0]\n";
print STDERR "ARG2 (NCOMP)  = $ARGV[1]\n";
print STDERR "------------------------------------------------------\n";
print STDERR "\n";
if ( $numArgs != 2 ) {
    print STDERR "ERROR: command line args, ($numArgs/2)\n";
    print STDERR "USAGE: GET_spect_components.pl [DIR] [NCOMP]\n";
    print STDERR "------------------------------------------------------\n";
    print STDERR "\n";
    exit(1);
}

$DIR   = $ARGV[0];
$NCOMP = $ARGV[1];

# Read the beamwidth from paramter file
open( BEAM_INFO, "${DIR}/${beamParamters}" )
  || die "beam parameter file not found!";
$DELTAF = <BEAM_INFO>;
$DELTAF = substr( $DELTAF, 0, 10 );    ## removes the endl from the string
close(BEAM_INFO);

# print a summary of the options
print STDOUT "======================================================\n";
print STDOUT "SPEC_FILE = ${DIR}/${spectrum}\n";
print STDOUT "NCOMP     = ${NCOMP}\n";
print STDOUT "======================================================\n";
print STDOUT "CBEAM FILE= ${DIR}/${beamParamters}\n";
print STDOUT "DELTAF    = $DELTAF\n";
print STDOUT "======================================================\n";

# Load the spectra from file
open( DATA, "${DIR}/${spectrum}" ) || die "spectrum file not found!";
@spect = <DATA>;
close(DATA);

# Extract values from the data (from start end end freq) 
$i = 0;
foreach (@spect) {
    split;
    if ($_[0] > 0 && $_[0] >= $startFreq && $_[0] <= $endFreq) {
        $freq[$i]   = $_[0];
        $amp[$i]    = $_[1];
        $phase[$i]  = $_[4];
        $cPhase[$i] = $_[6];
        $i++;
    }
}
$nfreqs = $i;

# sort the values in order of amplitude
@samp   = reverse sort { $a <=> $b } @amp;

$start   = 0;
$oldAMP  = 0;
$oldFREQ = 0;
$DIFF    = 0;
$NPEAKS  = 0;

print STDOUT "\n";
print STDOUT "freq             (period)         amp*2        phase/pi   corrected phase/pi    i     (iPeak) \n";
print STDOUT "============================================================================================\n";
for ( $peak = $start ; $peak < $nfreqs ; $peak++ ) {
    for ( $i = 0 ; $i < $nfreqs ; $i++ ) {
        $pass = 0;
        if ( $samp[$peak] == $amp[$i] ) {
            for ( $test = 0 ; $test < $NPEAKS ; $test++ ) {
                $min = $PEAKS[$test][1] - $DELTAF;
                $max = $PEAKS[$test][1] + $DELTAF;
                if ( $freq[$i] > $min && $freq[$i] < $max ) {
                    $pass = 1;
                }
            }

            if ( $pass == 0 ) {
                $PEAKS[$NPEAKS][0] = $i;
                $PEAKS[$NPEAKS][1] = $freq[$i];
                $PEAKS[$NPEAKS][2] = $amp[$i];
                printf(
                    "%3.6f\t(%3.6f)\t%3.6f\t%3.6f\t%3.6f\t%3d\t%3d\n",
                    $freq[$i],
                    1 / $freq[$i],
                    $amp[$i] * 2,
                    $phase[$i]/$PI,
                    $cPhase[$i]/$PI,
                    $i,
                    $NPEAKS+1
                );
                $NPEAKS++;
                if ( $NPEAKS == $NCOMP ) { exit(); }
            }
            $oldAMP  = $amp[$i];
            $oldFREQ = $freq[$i];
            $start++;
        }
    }
}

